
# Load packages
  library(tidyverse)
  library(car)
  library(lmtest)
  library(sandwich)
  library(stargazer)
  
# Calculate robust confidence intervals
  se_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Std. Error"]
  p_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Pr(>|t|)"]

  #### AUTNES Online Panel Study 2017 - 2019
  
  #### 2017
  
  ### wave 4
  load("ObsStudy_00_data_Austria_AUTNES_Online_Panel_2017_2019.RData")
  dt<- table

  # Party rating: Generate party rating as distance measure retrieved from left-right party placement and self-placement 
  dt[, 129:134][dt[, 129:134] >= 77] <- NA
  dt[, 135][dt[, 135] >= 77] <- NA
  
  # Take care of the placement variables for individuals that joined panel in wave 4
  dt[, 1042:1047][dt[, 1042:1047] >= 77] <- NA
  dt[, 1048][dt[, 1048] >= 77] <- NA
  
  # SPO, add information for individuals that joined panel in wave 4
  dt$w1_q8x1 <- ifelse(is.na(dt$w1_q8x1), dt$w4f_q7x1, dt$w1_q8x1)
  
  # OVP, add information for individuals that joined panel in wave 4
  dt$w1_q8x2 <- ifelse(is.na(dt$w1_q8x2), dt$w4f_q7x2, dt$w1_q8x2)
  
  # FPO, add information for individuals that joined panel in wave 4
  dt$w1_q8x3 <- ifelse(is.na(dt$w1_q8x3), dt$w4f_q7x3, dt$w1_q8x3)
  
  # GRUENE, add information for individuals that joined panel in wave 4
  dt$w1_q8x4 <- ifelse(is.na(dt$w1_q8x4), dt$w4f_q7x4, dt$w1_q8x4)
  
  # NEOS, add information for individuals that joined panel in wave 4
  dt$w1_q8x5 <- ifelse(is.na(dt$w1_q8x5), dt$w4f_q7x5, dt$w1_q8x5)
  
  # Stronach, add information for individuals that joined panel in wave 4
  dt$w1_q8x6 <- ifelse(is.na(dt$w1_q8x6), dt$w4f_q7x6, dt$w1_q8x6)
  
  # Self-placement, add information for individuals that joined panel in wave 4
  dt$w1_q9 <- ifelse(is.na(dt$w1_q9), dt$w4f_q8, dt$w1_q9)
  
  
  # Party rating SPO
  dt$rat_SPO <- abs(dt$w1_q9 - dt$w1_q8x1)^2
  # Recode, so large distance value is a low value on party rating variable
  dt$rat_SPO = recode(dt$rat_SPO, "0=10; 1=9; 2=8; 3=7; 4=6; 5=5; 6=4; 7=3; 8=2; 9=1; 10=0")

  # Party rating OVP
  dt$rat_OVP <- abs(dt$w1_q9 - dt$w1_q8x2)^2
  # Recode, so large distance value is a low value on party rating variable
  dt$rat_OVP = recode(dt$rat_OVP, "0=10; 1=9; 2=8; 3=7; 4=6; 5=5; 6=4; 7=3; 8=2; 9=1; 10=0")

  # Party rating FPO
  dt$rat_FPO <- abs(dt$w1_q9 - dt$w1_q8x3)^2
  # Recode, so large distance value is a low value on party rating variable
  dt$rat_OVP = recode(dt$rat_OVP, "0=10; 1=9; 2=8; 3=7; 4=6; 5=5; 6=4; 7=3; 8=2; 9=1; 10=0")

  # Party rating GRUENE
  dt$rat_GRUENE <- abs(dt$w1_q9 - dt$w1_q8x4)^2
  # Recode, so large distance value is a low value on party rating variable
  dt$rat_GRUENE = recode(dt$rat_GRUENE, "0=10; 1=9; 2=8; 3=7; 4=6; 5=5; 6=4; 7=3; 8=2; 9=1; 10=0")

  # Party rating NEOS
  dt$rat_NEOS <- abs(dt$w1_q9 - dt$w1_q8x5)^2
  # Recode, so large distance value is a low value on party rating variable
  dt$rat_NEOS = recode(dt$rat_NEOS, "0=10; 1=9; 2=8; 3=7; 4=6; 5=5; 6=4; 7=3; 8=2; 9=1; 10=0")
  
  # Coalition evaluations/preferred coalition
  # SPOE-OEVP, w4_q41x1 / OEVP-FPOE, w4_q41x2 / SPOE-FPOE, w4_q41x3 / OEVP-GRUENE-NEOS, w4_q42x1 / SPOE-GRUENE-NEOS, w4_q42x2 
  colnames(dt)[seq(824, 828, 1)] <- c("rat_coal_SPO_OVP", "rat_coal_OVP_FPO", "rat_coal_SPO_FPO", "rat_coal_OVP_GRUENE_NEOS", "rat_coal_SPO_GRUENE_NEOS")
  dt[, 824:828][dt[, 824:828] >= 77] <- NA
  
  # Coalition likelihood
  # SPOE-OEVP, w4_q40x1/ OEVP-FPOE, w4_q40x2 / SPOE-FPOE, w4_q40x3 
  colnames(dt)[seq(821, 823, 1)] <- c("coallik_SPO_OVP", "coallik_OVP_FPO", "coallik_SPO_FPO")
  dt[, 821:823][dt[, 821:823] >= 77] <- NA
  
  # Dependent Variable: binary vote choice
  #Q16.Und welcher Partei werden Sie bei der Nationalratswahl am 15. Oktober 2017 voraussichtlich Ihre Stimme geben?
  #Q18.Angenommen Sie w?rden an der Nationalratswahl am 15. Oktober 2017 teilnehmen, f?r welche Partei w?rden Sie sich am ehesten entscheiden?
  #Q19.Und welche Partei haben Sie bei der Nationalratswahl 2017 gew?hlt?
  
  dt$w4_q16[dt$w4_q16 >= 77] <- NA
  dt$w4_q18[dt$w4_q18 >= 77] <- NA
  dt$w4_q19[dt$w4_q19 >= 77] <- NA
  
  # Generate party choice variable
  dt$vote <- dt$w4_q16
  dt$vote <- ifelse(is.na(dt$vote), dt$w4_q18, dt$vote)
  dt$vote <- ifelse(is.na(dt$vote), dt$w4_q19, dt$vote)
  
  # Generate party choice columns based on vote variable
  dt$vc_SPO <- ifelse(dt$vote == 1, 1, 0)
  dt$vc_OVP <- ifelse(dt$vote == 2, 1, 0)
  dt$vc_FPO <- ifelse(dt$vote == 3, 1, 0)
  
  # Coalition evaluation
  Z <- dt %>% select(contains("rat_coal")) %>%
    mutate_all(as.numeric) %>%
    mutate_all(funs(scales::rescale(.,to = c(0, 1))))
  
  rat_coal_SPO <- Z %>% select(rat_coal_SPO_OVP, rat_coal_SPO_FPO) %>% as.matrix()
  rat_coal_OVP <- Z %>% select(rat_coal_SPO_OVP, rat_coal_OVP_FPO) %>% as.matrix()
  rat_coal_FPO <- Z %>% select(rat_coal_OVP_FPO, rat_coal_SPO_FPO) %>% as.matrix()
  
  # Coalition likelihood
  gamma <- dt %>% select(contains("coallik_")) 
  
  coal_lik_SPO <- gamma %>% select(coallik_SPO_OVP, coallik_SPO_FPO) %>%
    mutate(sum_exp = exp(coallik_SPO_OVP) + exp(coallik_SPO_FPO),
           coallik_SPO_OVP = exp(coallik_SPO_OVP)/sum_exp,
           coallik_SPO_FPO = exp(coallik_SPO_FPO)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  
  coal_lik_OVP <- gamma %>% select(coallik_SPO_OVP, coallik_OVP_FPO) %>%
    mutate(sum_exp = exp(coallik_SPO_OVP) + exp(coallik_OVP_FPO),
           coallik_SPO_OVP = exp(coallik_SPO_OVP) / sum_exp,
           coallik_OVP_FPO = exp(coallik_OVP_FPO) / sum_exp) %>%
    select(-sum_exp)  %>%
    as.matrix()
  
  
  coal_lik_FPO <- gamma %>% select(coallik_OVP_FPO, coallik_SPO_FPO) %>%
    mutate(sum_exp = exp(coallik_OVP_FPO) + exp(coallik_SPO_FPO),
           coallik_OVP_FPO = exp(coallik_OVP_FPO) / sum_exp,
           coallik_SPO_FPO = exp(coallik_SPO_FPO) / sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  
  # Mean-variance model
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  N <- nrow(dt)
  
  M_SPO <- sapply(1:N, function(i) EV(coal_lik_SPO[i,], rat_coal_SPO[i,]))
  V_SPO <- sapply(1:N, function(i) Vari(coal_lik_SPO[i,], rat_coal_SPO[i,]))
  
  M_OVP <- sapply(1:N, function(i) EV(coal_lik_OVP[i,], rat_coal_OVP[i,]))
  V_OVP <- sapply(1:N, function(i) Vari(coal_lik_OVP[i,], rat_coal_OVP[i,]))
  
  M_FPO <- sapply(1:N, function(i) EV(coal_lik_FPO[i,], rat_coal_FPO[i,]))
  V_FPO <- sapply(1:N, function(i) Vari(coal_lik_FPO[i,], rat_coal_FPO[i,]))
  
  
  # Gender variable, 1 = male, 2 = female
  dt$male <- ifelse(dt$sd3 == 2, 0, dt$sd3)
  dt$male <- ifelse(dt$male==3, NA, dt$male)

  # Education
  # Set missings or observations that do not fit the scheme to NA
  dt$sd7 <- ifelse(dt$sd7 > 15, NA, dt$sd7)
  
  # PID
  if (!exists("robustnesscheck_pid")) {
    robustnesscheck_pid <- FALSE # PID as robustness?
  }
  dt$w1_pid <- dt$w1_q19
  dt$w1_pid[dt$w1_q18==2] <- 0
  dt$w4_pid <- dt$w4f_q26
  dt$w4_pid[dt$w4f_q25==2] <- 0
  dt$pid <- ifelse(!is.na(dt$w4_pid), dt$w4_pid, dt$w1_pid)
  
  dt$pid_SPO <- ifelse(dt$pid==1, 1, 0)
  dt$pid_OVP <- ifelse(dt$pid==2, 1, 0)
  dt$pid_FPO <- ifelse(dt$pid==3, 1, 0)

  d <- data.frame(
    "vc_SPO" = dt$vc_SPO,
    "vc_OVP" = dt$vc_OVP,
    "vc_FPO" = dt$vc_FPO,
    "rat.SPO" = dt$rat_SPO,
    "rat.OVP" = dt$rat_OVP,
    "rat.FPO" = dt$rat_FPO,
    "pid_SPO" = dt$pid_SPO,
    "pid_OVP" = dt$pid_OVP,
    "pid_FPO" = dt$pid_FPO,
    "lotterymean.SPO" = M_SPO,
    "lotteryvariance.SPO" = V_SPO,
    "lotterymean.OVP" = M_OVP,
    "lotteryvariance.OVP" = V_OVP,
    "lotterymean.FPO" = M_FPO,
    "lotteryvariance.FPO" = V_FPO,
    "sex" = dt$male,
    "age" = dt$age,
    "edu" = dt$sd7,
    "rat_leader_SPO" = ifelse(dt$w4_q11x1<11, dt$w4_q11x1, NA),
    "rat_leader_OVP" = ifelse(dt$w4_q11x2<11, dt$w4_q11x2, NA),
    "rat_leader_FPO" = ifelse(dt$w4_q11x3<11, dt$w4_q11x3, NA)
  )  
  
  d$lotteryvariance.SPO <- scales::rescale(d$lotteryvariance.SPO, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.OVP <- scales::rescale(d$lotteryvariance.OVP, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.FPO <- scales::rescale(d$lotteryvariance.FPO, to = c(0, 1), from = c(0,0.25))
  
# Save model results
  summary(m1 <- lm(as.formula(paste("vc_SPO ~ lotteryvariance.SPO + lotterymean.SPO + rat.SPO + rat_leader_SPO + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_SPO" else "")), d))
  summary(m2 <- lm(as.formula(paste("vc_OVP ~ lotteryvariance.OVP + lotterymean.OVP + rat.OVP + rat_leader_OVP + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_OVP" else "")), d))
  summary(m3 <- lm(as.formula(paste("vc_FPO ~ lotteryvariance.FPO + lotterymean.FPO + rat.FPO + rat_leader_FPO + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_FPO" else "")), d))
  
  AUT_2017_AUTNES_OP_w4_SPO_Variance_Estimate <- tidy(m1) %>% filter(term=="lotteryvariance.SPO") %>% select("estimate")
  AUT_2017_AUTNES_OP_w4_SPO_Variance_SE <- se_robust(m1)["lotteryvariance.SPO"]
  AUT_2017_AUTNES_OP_w4_SPO_Mean_Estimate <- tidy(m1) %>% filter(term=="lotterymean.SPO") %>% select("estimate")
  AUT_2017_AUTNES_OP_w4_SPO_Mean_SE <- se_robust(m1)["lotterymean.SPO"]
  
  AUT_2017_AUTNES_OP_w4_OVP_Variance_Estimate <- tidy(m2) %>% filter(term=="lotteryvariance.OVP") %>% select("estimate")
  AUT_2017_AUTNES_OP_w4_OVP_Variance_SE <- se_robust(m2)["lotteryvariance.OVP"]
  AUT_2017_AUTNES_OP_w4_OVP_Mean_Estimate <- tidy(m2) %>% filter(term=="lotterymean.OVP") %>% select("estimate")
  AUT_2017_AUTNES_OP_w4_OVP_Mean_SE <- se_robust(m2)["lotterymean.OVP"]
  
  AUT_2017_AUTNES_OP_w4_FPO_Variance_Estimate <- tidy(m3) %>% filter(term=="lotteryvariance.FPO") %>% select("estimate")
  AUT_2017_AUTNES_OP_w4_FPO_Variance_SE <- se_robust(m3)["lotteryvariance.FPO"]
  AUT_2017_AUTNES_OP_w4_FPO_Mean_Estimate <- tidy(m3) %>% filter(term=="lotterymean.FPO") %>% select("estimate")
  AUT_2017_AUTNES_OP_w4_FPO_Mean_SE <- se_robust(m3)["lotterymean.FPO"]
  
  # Harmonize names of the IVs of interest in the models
  names(m1$coefficients)[names(m1$coefficients) == "lotteryvariance.SPO"] <- "Government Lottery Variance"
  names(m1$coefficients)[names(m1$coefficients) == "lotterymean.SPO"] <- "Government Lottery Mean"
  names(m2$coefficients)[names(m2$coefficients) == "lotteryvariance.OVP"] <- "Government Lottery Variance"
  names(m2$coefficients)[names(m2$coefficients) == "lotterymean.OVP"] <- "Government Lottery Mean"
  names(m3$coefficients)[names(m3$coefficients) == "lotteryvariance.FPO"] <- "Government Lottery Variance"
  names(m3$coefficients)[names(m3$coefficients) == "lotterymean.FPO"] <- "Government Lottery Mean"
  
  
  if(!robustnesscheck_pid){
  stargazer(m1, m2, m3, title="Linear regressions of vote choice on perceived government lottery variance and mean (AUTNES Online Panel Study 2017-2019, wave 4, 2017)", 
            align=TRUE, 
            dep.var.labels=c("SPO", "OVP", "FPO"),
            omit.stat=c("LL","ser","f"),
            se = lapply(list(m1, m2, m3), se_robust),
            p = lapply(list(m1, m2, m3), p_robust),
            out = "TableSM4.tex")
  }
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  ### wave 11, pre-election 2019
  load("ObsStudy_00_data_Austria_AUTNES_Online_Panel_2017_2019.RData")
  dt<- table

  # Party rating: Generate party rating as distance measure retrieved from left-right party placement and self-placement 
  dt[, 129:134][dt[, 129:134] >= 77] <- NA
  dt[, 135][dt[, 135] >= 77] <- NA
  
  # Take care of the placement variables for individuals that joined panel in wave 11
  dt[, 2540:2544][dt[, 2540:2544] >= 77] <- NA
  dt[, 2546][dt[, 2546] >= 77] <- NA
  
  
  # SPO, add information for individuals that joined panel in wave 11
  dt$w1_q8x1 <- ifelse(is.na(dt$w1_q8x1), dt$w11f_q6x1, dt$w1_q8x1)
  
  # OVP, add information for individuals that joined panel in wave 11
  dt$w1_q8x2 <- ifelse(is.na(dt$w1_q8x2), dt$w11f_q6x2, dt$w1_q8x2)
  
  # FPO, add information for individuals that joined panel in wave 11
  dt$w1_q8x3 <- ifelse(is.na(dt$w1_q8x3), dt$w11f_q6x3, dt$w1_q8x3)
  
  # GRUENE, add information for individuals that joined panel in wave 11
  dt$w1_q8x4 <- ifelse(is.na(dt$w1_q8x4), dt$w11f_q6x4, dt$w1_q8x4)
  
  # NEOS, add information for individuals that joined panel in wave 11
  dt$w1_q8x5 <- ifelse(is.na(dt$w1_q8x5), dt$w11f_q6x5, dt$w1_q8x5)
  
  # Self-placement, add information for individuals that joined panel in wave 11
  dt$w1_q9 <- ifelse(is.na(dt$w1_q9), dt$w11f_q7, dt$w1_q9)
  
  
  # Party rating SPO
  dt$rat_SPO <- abs(dt$w1_q9 - dt$w1_q8x1)^2
  # Recode, so large distance value is a low value on party rating variable
  dt$rat_SPO = recode(dt$rat_SPO, "0=10; 1=9; 2=8; 3=7; 4=6; 5=5; 6=4; 7=3; 8=2; 9=1; 10=0")

  # Party rating OVP
  dt$rat_OVP <- abs(dt$w1_q9 - dt$w1_q8x2)^2
  # Recode, so large distance value is a low value on party rating variable
  dt$rat_OVP = recode(dt$rat_OVP, "0=10; 1=9; 2=8; 3=7; 4=6; 5=5; 6=4; 7=3; 8=2; 9=1; 10=0")

  # Party rating FPO
  dt$rat_FPO <- abs(dt$w1_q9 - dt$w1_q8x3)^2
  # Recode, so large distance value is a low value on party rating variable
  dt$rat_FPO = recode(dt$rat_FPO, "0=10; 1=9; 2=8; 3=7; 4=6; 5=5; 6=4; 7=3; 8=2; 9=1; 10=0")
  
  # Coalition evaluations/preferred coalition
  # SPOE-OEVP, w11_q36x1 / OEVP-FPOE, w11_q36x2 / SPOE-FPOE, w11_q36x3 / OEVP-NEOS, w11_q36x4 / OEVP-GRUENE-NEOS, w11_q37x1 / SPOE-GRUENE-NEOS, w11_q37x2
  colnames(dt)[seq(2662, 2665, 1)] <- c("rat_coal_SPO_OVP", "rat_coal_OVP_FPO", "rat_coal_SPO_FPO", "rat_coal_OVP_NEOS") # Don´t include OEVP minority government
  colnames(dt)[seq(2667, 2668, 1)] <- c("rat_coal_OVP_GRUENE_NEOS", "rat_coal_SPOE_GRUENE_NEOS")
  dt[, 2662:2668][dt[, 2662:2668] >= 77] <- NA
  
  # Coalition likelihood
  # SPOE-OEVP, w11_q35x1/ OEVP-FPOE, w11_q35x2 / SPOE-FPOE, w11_q35x3 / OEVP-GRUENE-NEOS, w11_q35x4 / SPOE-GRUENE-NEOS, w11_q35x5
  colnames(dt)[seq(2655, 2660, 1)] <- c("coallik_SPO_OVP", "coallik_OVP_FPO", "coallik_SPO_FPO", "coallik_OEVP_GRUENE_NEOS", "coallik_SPOE_GRUENE_NEOS", "coallik_OEVP_NEOS")
  dt[, 2655:2660][dt[, 2655:2660] >= 77] <- NA
  
  # Dependent variable: binary vote choice
  #Q15.Und welcher Partei werden Sie bei der Nationalratswahl am 29. September 2019 voraussichtlich Ihre Stimme geben?
  #Q17.Angenommen Sie w?rden an der Nationalratswahl am 29. September 2019 teilnehmen, f?r welche Partei w?rden Sie sich am ehesten entscheiden?
  #Q18.Und welche Partei haben Sie bei der Nationalratswahl 2019 gew?hlt?
  
  dt$w11_q15[dt$w11_q15 >= 77] <- NA
  dt$w11_q17[dt$w11_q17 >= 77] <- NA
  dt$w11_q18[dt$w11_q18 >= 77] <- NA
  
  # Generate party choice variable
  dt$vote <- dt$w11_q15
  dt$vote <- ifelse(is.na(dt$vote), dt$w11_q17, dt$vote)
  dt$vote <- ifelse(is.na(dt$vote), dt$w11_q18, dt$vote)
  
  # Generate party choice columns based on vote variable
  dt$vc_OVP <- ifelse(dt$vote == 1, 1, 0)
  dt$vc_SPO <- ifelse(dt$vote == 2, 1, 0)
  dt$vc_FPO <- ifelse(dt$vote == 3, 1, 0)
  
    
  # Coalition evaluation
  Z <- dt %>% select(contains("rat_coal")) %>%
    mutate_all(as.numeric) %>%
    mutate_all(funs(scales::rescale(.,to = c(0, 1))))
  
  rat_coal_SPO <- Z %>% select(rat_coal_SPO_OVP, rat_coal_SPO_FPO) %>% as.matrix()
  rat_coal_OVP <- Z %>% select(rat_coal_SPO_OVP, rat_coal_OVP_FPO) %>% as.matrix()
  rat_coal_FPO <- Z %>% select(rat_coal_OVP_FPO, rat_coal_SPO_FPO) %>% as.matrix()
  
  # Coalition likelihood
  gamma <- dt %>% select(contains("coallik_")) 
  
  coal_lik_SPO <- gamma %>% select(coallik_SPO_OVP, coallik_SPO_FPO) %>%
    mutate(sum_exp = exp(coallik_SPO_OVP) + exp(coallik_SPO_FPO),
           coallik_SPO_OVP = exp(coallik_SPO_OVP)/sum_exp,
           coallik_SPO_FPO = exp(coallik_SPO_FPO)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  
  coal_lik_OVP <- gamma %>% select(coallik_SPO_OVP, coallik_OVP_FPO) %>%
    mutate(sum_exp = exp(coallik_SPO_OVP) + exp(coallik_OVP_FPO),
           coallik_SPO_OVP = exp(coallik_SPO_OVP) / sum_exp,
           coallik_OVP_FPO = exp(coallik_OVP_FPO) / sum_exp) %>%
    select(-sum_exp)  %>%
    as.matrix()
  
  
  coal_lik_FPO <- gamma %>% select(coallik_OVP_FPO, coallik_SPO_FPO) %>%
    mutate(sum_exp = exp(coallik_OVP_FPO) + exp(coallik_SPO_FPO),
           coallik_OVP_FPO = exp(coallik_OVP_FPO) / sum_exp,
           coallik_SPO_FPO = exp(coallik_SPO_FPO) / sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  N <- nrow(dt)
  
  M_SPO <- sapply(1:N, function(i) EV(coal_lik_SPO[i,], rat_coal_SPO[i,]))
  V_SPO <- sapply(1:N, function(i) Vari(coal_lik_SPO[i,], rat_coal_SPO[i,]))
  
  M_OVP <- sapply(1:N, function(i) EV(coal_lik_OVP[i,], rat_coal_OVP[i,]))
  V_OVP <- sapply(1:N, function(i) Vari(coal_lik_OVP[i,], rat_coal_OVP[i,]))
  
  M_FPO <- sapply(1:N, function(i) EV(coal_lik_FPO[i,], rat_coal_FPO[i,]))
  V_FPO <- sapply(1:N, function(i) Vari(coal_lik_FPO[i,], rat_coal_FPO[i,]))
  
  
  # Gender variable, 1 = male, 2 = female
  dt$male <- ifelse(dt$sd3 == 2, 0, dt$sd3)
  dt$male <- ifelse(dt$male== 3, NA, dt$male)

  # Education
  # Set missings or observations that do not fit the scheme to NA
  dt$sd7 <- ifelse(dt$sd7 > 15, NA, dt$sd7)
  
  # PID
  if (!exists("robustnesscheck_pid")) {
    robustnesscheck_pid <- FALSE # PID as robustness?
  }
  dt$w7_pid <- dt$w7_q44
  dt$w7_pid[dt$w7_q43==2] <- 0
  
  dt$w8_pid <- dt$w8_q18
  dt$w8_pid[dt$w8_q17==2] <- 0
  
  dt$w9_pid <- dt$w9f_q57
  dt$w9_pid[dt$w9f_q56==2] <- 0
  
  dt$w10_pid <- dt$w10f_q84
  dt$w10_pid[dt$w10f_q83==2] <- 0
  
  dt$w11_pid <- dt$w11f_q25
  dt$w11_pid[dt$w11f_q24==2] <- 0
  
  dt$pid <- ifelse(!is.na(dt$w11_pid), dt$w11_pid, dt$w10_pid)
  dt$pid <- ifelse(!is.na(dt$pid), dt$pid, dt$w9_pid)
  dt$pid <- ifelse(!is.na(dt$pid), dt$pid, dt$w8_pid)
  dt$pid <- ifelse(!is.na(dt$pid), dt$pid, dt$w7_pid)
  
  dt$pid_SPO <- ifelse(dt$pid==2, 1, 0)
  dt$pid_OVP <- ifelse(dt$pid==1, 1, 0)
  dt$pid_FPO <- ifelse(dt$pid==3, 1, 0)
  
  d <- data.frame(
    "vc_SPO" = dt$vc_SPO,
    "vc_OVP" = dt$vc_OVP,
    "vc_FPO" = dt$vc_FPO,
    "rat.SPO" = dt$rat_SPO,
    "rat.OVP" = dt$rat_OVP,
    "rat.FPO" = dt$rat_FPO,
    "pid_SPO" = dt$pid_SPO,
    "pid_OVP" = dt$pid_OVP,
    "pid_FPO" = dt$pid_FPO,
    "lotterymean.SPO" = M_SPO,
    "lotteryvariance.SPO" = V_SPO,
    "lotterymean.OVP" = M_OVP,
    "lotteryvariance.OVP" = V_OVP,
    "lotterymean.FPO" = M_FPO,
    "lotteryvariance.FPO" = V_FPO,
    "sex" = dt$male,
    "age" = dt$age,
    "edu" = dt$sd7,
    "rat_leader_SPO" = ifelse(dt$w11_q11x1<11, dt$w11_q11x1, NA),
    "rat_leader_OVP" = ifelse(dt$w11_q11x2<11, dt$w11_q11x2, NA),
    "rat_leader_FPO" = ifelse(dt$w11_q11x3<11, dt$w11_q11x3, NA)
  )  
  
  d$lotteryvariance.SPO <- scales::rescale(d$lotteryvariance.SPO, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.OVP <- scales::rescale(d$lotteryvariance.OVP, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.FPO <- scales::rescale(d$lotteryvariance.FPO, to = c(0, 1), from = c(0,0.25))
  
# Save model results
  summary(m1 <- lm(as.formula(paste("vc_SPO ~ lotteryvariance.SPO + lotterymean.SPO + rat.SPO + rat_leader_SPO + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_SPO" else "")), d))
  summary(m2 <- lm(as.formula(paste("vc_OVP ~ lotteryvariance.OVP + lotterymean.OVP + rat.OVP + rat_leader_OVP + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_OVP" else "")), d))
  summary(m3 <- lm(as.formula(paste("vc_FPO ~ lotteryvariance.FPO + lotterymean.FPO + rat.FPO + rat_leader_FPO + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_FPO" else "")), d))
  
  AUT_2019_AUTNES_OP_w11_SPO_Variance_Estimate <- tidy(m1) %>% filter(term=="lotteryvariance.SPO") %>% select("estimate")
  AUT_2019_AUTNES_OP_w11_SPO_Variance_SE <- se_robust(m1)["lotteryvariance.SPO"]
  AUT_2019_AUTNES_OP_w11_SPO_Mean_Estimate <- tidy(m1) %>% filter(term=="lotterymean.SPO") %>% select("estimate")
  AUT_2019_AUTNES_OP_w11_SPO_Mean_SE <- se_robust(m1)["lotterymean.SPO"]
  
  AUT_2019_AUTNES_OP_w11_OVP_Variance_Estimate <- tidy(m2) %>% filter(term=="lotteryvariance.OVP") %>% select("estimate")
  AUT_2019_AUTNES_OP_w11_OVP_Variance_SE <- se_robust(m2)["lotteryvariance.OVP"]
  AUT_2019_AUTNES_OP_w11_OVP_Mean_Estimate <- tidy(m2) %>% filter(term=="lotterymean.OVP") %>% select("estimate")
  AUT_2019_AUTNES_OP_w11_OVP_Mean_SE <- se_robust(m2)["lotterymean.OVP"]
  
  AUT_2019_AUTNES_OP_w11_FPO_Variance_Estimate <- tidy(m3) %>% filter(term=="lotteryvariance.FPO") %>% select("estimate")
  AUT_2019_AUTNES_OP_w11_FPO_Variance_SE <- se_robust(m3)["lotteryvariance.FPO"]
  AUT_2019_AUTNES_OP_w11_FPO_Mean_Estimate <- tidy(m3) %>% filter(term=="lotterymean.FPO") %>% select("estimate")
  AUT_2019_AUTNES_OP_w11_FPO_Mean_SE <- se_robust(m3)["lotterymean.FPO"]
  
  # Harmonize names of the IVs of interest in the models
  names(m1$coefficients)[names(m1$coefficients) == "lotteryvariance.SPO"] <- "Government Lottery Variance"
  names(m1$coefficients)[names(m1$coefficients) == "lotterymean.SPO"] <- "Government Lottery Mean"
  names(m2$coefficients)[names(m2$coefficients) == "lotteryvariance.OVP"] <- "Government Lottery Variance"
  names(m2$coefficients)[names(m2$coefficients) == "lotterymean.OVP"] <- "Government Lottery Mean"
  names(m3$coefficients)[names(m3$coefficients) == "lotteryvariance.FPO"] <- "Government Lottery Variance"
  names(m3$coefficients)[names(m3$coefficients) == "lotterymean.FPO"] <- "Government Lottery Mean"
  
  if(!robustnesscheck_pid){
  stargazer(m1, m2, m3, title="Linear regressions of vote choice on perceived government lottery variance and mean (AUTNES Online Panel Study 2017-2019, wave 11, 2019)", 
            align=TRUE, 
            dep.var.labels=c("SPO", "OVP", "FPO", "GRUENE", "NEOS"),
            omit.stat=c("LL","ser","f"),
            type = "text",
            out = "TableSM5.tex"
            )
}